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ABSTRACT 

We develop a computational method for the estimation of parameters in a 
distributed model for a flexible structure. The structure we consider (part 
of the "RPL experiment") consists of a cantilevered beam with a thruster and 
linear accelerometer at the free end. The thruster is fed by a pressurized 
hose whose horizontal motion effects the transverse vibration of the beam. We 
use the Euler-Bernoulli theory to model the vibration of the beam and treat 
the hose-thruster assembly as a lumped or point mass-dashpot-spring system at 
the tip. Using measurements of linear acceleration at the tip, we estimate 
the hose parameters (mass, stiffness, damping) and a Voigt-Kelvin viscoelastic 
structural damping parameter for the beam using a least squares fit to the 
data . 


We consider spline based approximations to the hybrid (coupled ordinary 
and partial differential equations) system; theoretical convergence results 
and numerical studies with both simulation and actual experimental data 
obtained from the structure are presented and discussed. 
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1. Introduction 

The difficulties involved in the design of practical and effi- 
cient control laws for large flexible spacecraft (e.g. the inherent 
infinite dimensionality of the system, a large number of closely 
spaced modal frequencies, high flexibility, light damping, a fuel- 
limited , hostile , highly variable environment , etc . ) have stimulated 
research into the development of system identification and parameter 
estimation procedures which will yield high fidelity models. A partic- 
ular area of interest involves schemes for the estimation of material 
parameters describing, for example, mass, inertia, stiffness or 
damping properties in distributed models for the vibration of 
viscoelastic systems-specifically , mechanical beams, plates and the 
like. In addition, since the resulting inverse problems are often 
infinite dimensional, substantial attention has been focused on 
approximation; see, for example, [1], [2], [3], [4], [8] and [12]. 

In these treatments, the parameter estimation problem is formulated as 
a least squares fit to measurements of either displacement or 
velocity. Although significant gains have been made in the development 
of instrumentation to measure displacement and velocity (e.g. laser 
technology, etc.), one of the least expensive, most reliable and most 
commonly used sensors is the linear accelerometer . While in principle 
it is possible to integrate acceleration measurements once or twice 
to obtain respectively velocity or displacement data, in practice this 
task can pose significant challenges. For example, integration of the 
signal could result in the amplification of low frequency measurement 
noise or dynamic effects which have not been included in the underly- 
ing model. In light of this, we have undertaken to show here, both 
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theoretically and computationally, that a scheme in the spirit of 
those developed in the previously cited references can also be 
effectively used with acceleration measurements. In particular we 
note, this involves the nontrivial extension of the familiar 
variational arguments which are used to demonstrate the convergence of 
the finite element state approximations upon which the identification 
schemes are based. Indeed, it must be shown that in addition to the 
convergence of the displacement and velocity, the convergence of 
acceleration can be obtained as well. 

The other primary motivation for the present effort is that while 
these methods have been extensively tested and evaluated with simula- 
tion data, they have never been tried with actual experimental data. 

We have tested our scheme with data obtained from an experimental 
structure which was designed and constructed at the Charles Stark 
Draper Laboratory in Cambridge, Massachusetts with funding provided by 
the United States Air Force Rocket Propulsion Laboratory (RPL) . The 
RPL structure (as it will henceforth be referred to as) was designed 
to serve as a test bed for the implementation and evaluation of 
control algorithms for large angle slewing of spacecraft with flexible 
appendages. The structure was specifically designed to exhibit 
structural modes and damping characteristics representative of 
realistic large flexible space structures . 

In Section 2 we describe the RPL structure (its geometry, 
instrumentation, etc.) and formulate an inverse problem involving a 
distributed system. In Section 3, we use the resulting infinite 
dimensional estimation problem to motivate the development of a finite 
dimensional, finite element based approximation scheme. We also 
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discuss our theoretical convergence results. In Section 4 we present 
numerical findings. 

We use standard notation throughout. For X a normed linear 
space, L(X) denotes the space of bounded linear operators from X into 

V 

X. For Cl an interval and k = 0,1,2,--*, C (fl;X) denotes the space of 
functions from n into X which are k times continuously strongly 
differentiable on Cl. When k = 0 we shall simply write C(n;X). A 
function f from Cl into X will be said to belong to L 2 (fl;X) if 
Ifi|f (t) |x * 00 ' For k = 0,1,2, •••, H k (fl;X) denotes the completion 

of C (0;X) with respect to the norm 


Mk - 'j 0 In 

If, in addition, X is a Hilbert space with inner product < • , • > x then 
H (fi;X)is a Hilbert space with inner product 




l f n <f ( ^(t), g (J) (t)>^ dt. 
J-0 Jn X 


Tr 

When X = R, we use the abbreviated notations C (n), L 2 (n) and H (n) . 
Note that H # (n) = L 2 (0) and <*,*>q is the standard inner product on 

L,(n). 
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2. The Identification Problem 

The RPL structure (see Figure 2.1 below) consists of four 
flexible appendages which are cantilevered at right angles to one 
another from a rigid central hub. The hub is mounted on an air 
bearing table thus permitting the near frictionless rotation of the 
structure about the vertical axis. 


FLEXIBLE HOSES 



Figure 2 . 1 

Two of the appendages (which are mounted to the hub 180° apart) are 
"active"; each has two nitrogen cold gas thrusters mounted in opposing 
directions at its tip. The remaining two appendages are "passive" 
with only counter-balancing masses affixed to their free ends . The 
presence of the tip masses on the passive arms serves to preserve the 
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overall symmetry of the structure. Nitrogen gas from tanks mounted on 
the central hub is supplied to the thrusters via two stainless steel 
mesh-wrapped high pressure hoses. The expulsion of propellant from 
the thruster nozzles is controlled by electro-mechanical or solenoidal 
valves. Each of the four appendages is equipped with a sensor in the 
form of a linear accelerometer attached at its tip. Data from the 
accelerometers is processed and recorded and control input signals to 
the thrusters are generated by a MINC 11/23 microcomputer. A detailed 
description of the structure's design specifications can be found in 
[ 6 ] and [15]. 

The problem which is of primary concern to us here involves the 
modeling of the effects of the nitrogen supply hoses on the transverse 
vibration of the active members. We consider therefore, the structure 
with the central hub immobilized and look only at the vibration of one 
of the active appendages and view it as a simple cantilevered beam 
(see Figure 2.2). 



Figure 2.2 
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We treat the thruster assembly as a point mass that is rigidly attach- 
ed to the beam at the tip and propose a model for the hose effects in 
the form of a proof mass which reacts against the tip mass . In 
effect, we consider the idealized, simplified structure depicted in 
Figure 2.3 below involving a single, cantilevered, flexible, uniform 
beam with a two-mass-dashpot-spring system affixed to its free end. 





Figure 2.3 

In formulating a mathematical model for the structure shown in 
Figure 2.3 above, we assume that the beam is of length J. with uniform 
rectangular cross section of height h and width b. We let u(t,x) and 
y(t) denote respectively the transverse displacement of the beam at 
position x along its span and the displacement of the proof or hose 
mass, each at time t. Both are measured relative to the x-axis in the 
coordinate frame determined by the longitudinal axis of the beam in 
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its undeformed state with origin located at the beam's root or fixed 
end. Assuming the beam undergoes only small deformations (i.e. 

|u(t,x)| << Jl and |hp@F — ( t,x)| << 1) and has a small height to span 

length ratio, the Euler-Bernoulli theory (see [5]) including Voigt- 
Kelvin viscoelastic structural damping (see [10]) yields the partial 
differential equation 

8*u 8* 8u 8 4 u 

(2.1) p (t,x) + o n I (t,x) + El (t.x) « 0, 

8t* 8x 4 8t 8x 4 

0 < x < i, t > 0 

where p is the linear mass density of the beam, E is the modulus of 
elasticity, c D is the coefficient of viscosity and I is the second 
moment or moment of inertia of the cross sectional area A about the 
neutral axis. For the beam we consider here with constant rectangular 
cross section, I - bh 3 /12. Since the beam is assumed to be uniform, 
the parameters p, E and c^ are taken to be constant in time and space. 

Balancing forces at the free end, elementary Newtonian mechanics 
yields the equations of motion 

( 2 . 2 ) 

and 
(2.3) 


8*u 8* 8u 8 3 u 

m_ — (t.JL) - C-.I — — (t,D - El — (t.X.) 

1 8t * u 8x 3 8t 8x 3 


dv 3u 

c w (— (t) (t.JL)) + (y(t) - u(t,JL)) + f(t), 

H dt 8t H 


t > 0 


m w ^-^(t) + c„ (— (t) - — (t.i)) + k„ (y(t) - u(t , i) ) = 0, 
H dt* H dt 9t 


t > 0 



for the tip and hose masses m T and m H respectively. Here k H is the 
hose stiffness, c H is the hose damping coefficient and f(t) is the 
externally applied force at time t due to the firing of the thrusters 
mounted at the tip. 

Making the assumption that the rotatory inertia of the proof mass 
system is negligible, rotational equilibrium at the tip can be 
expressed as 

3 2 9u 3 2 u 

(2.4) o_.I (t,l) + El (t , JL) =0, t > 0. 

3x 2 3t 3x 2 

The zero displacement and zero slope constraints at the fixed end are 
given by 


(2.5) u(t , 0) = 0 


and 


3u 

— (t,0) = 0, 
3x 


t > 0 


respectively. Taking the structure to be initially at rest we have 
the initial conditions 


(2.6) u(0,x) = 0 and 

and 


3u 

— (0,x) = 0, 
3t 


(2.7) y(0) = 0 


and 
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In the mathematical model given by (2.1) - (2.7) above the parameters 
p , m T and I can be measured or computed directly . The modulus of 
elasticity E is typically determined in the laboratory. For the most 
commonly used materials (including aluminum which is the material from 
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which the structure of interest to us here is made) its value can be 
readily looked up in standard engineering tables. The parameters c D> 
m H’ c h an< * k H on k an< ** roust be determined experimentally; 

that is, they will have to be identified based upon the observed 
response of the structure to a given input disturbance. This is one 
class of inverse problems which we formulate and consider below. In 
the system of equations (2.1) - (2.7) we explicitly modeled (albeit, 
in a rather simple fashion) the dynamical effects of the hose. The 
unknown hose parameters are then determined as the solution to an 
inverse problem. 

An alternative approach to obtaining a model which exhibits a 
reasonable degree of fidelity involves a technique whioh is sometimes 
referred to as model adjustment. Starting with a simple model, the 
parameters are then "adjusted" so as to compensate for unraodeled 
dynamics. The choice of parameters to be adjusted and the resulting 
variations may or may not be motivated by physical considerations. 

In our problem for example, we might consider a simple cantile- 
vered beam with tip mass (i.e. m^ «= c^ = k fi = 0) and then adjust the 
theoretical or measured values of E and m T to compensate for the 
dynamical effects which result from the hose mass and motion. A value 
for the parameter c D could also be identified if damping effects are 
considered significant. Model adjustment was used in [6] to obtain a 
model for the RPL structure upon which control design could be based. 

We define an inverse problem which encompasses both of the 
general approaches which have been outlined above. We assume that an 
input disturbance described by the function f(t), t e [0,T] is applied 
to the structure via the tip thrusters and that the linear accelera- 
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tion at the free end of the beam, z(t), is measured and recorded for 
each t e [t,,^] where O^to^ti^T. (Of course, in actual 
practice, z could in fact only be sampled discretely). Let R + denote 
the positive real numbers and let Q be a closed and bounded subset of 
R°. We seek a q e Q which minimizes 

J(q) “ |“ 7 (t,X,;q) - z(t)| 2 dt 

where u(*,*;q) denotes the solution to the initial-boundary value 
problem (2.1) - (2.7) corresponding to q = (m T ,E ,c D ,m H , c H ,k H ) e Q. 

Our primary concerns in the next section will include well- 
posedness of the system (2.1) - (2.7), existence of a minimizer for J, 
and development of approximation techniques to find this minimizer. 
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3. Approximation Theory 

A computational method for the solution of the estimation problem 
posed above will invariably involve finite dimensional approximation 
of the initial -boundary value problem (2.1) - (2.7). We have been 
successful in solving inverse problems for distributed parameter 
models for flexible structures (see, for example, [1], [2], [3], [4], 
[12]) using spline-based Ritz-Galerkin techniques. We apply those 
ideas here and derive finite element approximations based upon an 
abstract Hilbert space formulation of the hybrid system of ordinary 
and partial differential equations and boundary conditions given in 
(2.1) - (2.7). This abstract formulation is also useful in 
establishing existence, uniqueness and necessary regularity results 
for solutions. We briefly outline the essential features of our 
general approach (including theoretical convergence results) in the 
context of the particular problem of interest to us here. 

Let H - R* x L 2 (0 , X.) be endowed with the usual product space 
inner product 

<(S,H,<IO.(X,U,nO> h - ?X + mi + <$,V > 0 


and let 


v - (U,P,<D) € H: 4> € H* (0 , X.) , 4>(0) - D0(O) - 0, r\ « <D(X,)} 


be endowed with the inner product 


<(?,0U),0),(X,N>(X.),lO> v = (?-<t>(X.) ) (X-N'(X-) ) + <D 2 <t>,D 2 N>> 0 


12 


where the symbol D is used here and below to denote the spatial 

d 

differentiation operator — . The space V together with the inner 

dx 

product < • , • > v form a Hilbert space which is densely and compactly 
embedded in H. 

We rewrite the system (2.1) - (2.7) as the abstract second order 
initial value problem in H 

(3.1) Mu tt (t) + Cu t (t) + Ku(t) - F(t) , t > 0 

(3.2) 6(u(t) + eu t (t)) =0, t > 0 

(3.3) u(0) = 0 u t (0) * 0 

a. 

in the states u(t) = (y(t) ,u(t , i) ,u(t, • )) . The operators M € L(H) , 
C:D C H -+ H and K:D C H -» H are given by 

(3.4) M(^ , -n , <()) = (m H ^,m T n,p<t>), 

cu,n,<t») = (c H (^-Ti),c H ('n-;) - o D i D 3 <t>a),c D i d 4 4>), 

and 

KU.n.dO = (k H U-n),k H (n-0 - EI d’KjO.ei d 4 <i>) 

where D = {(^ ,r» , «t>) € V: <t> € H 4 (0,!l)}. For each t > 0, F(t) = 
(0,f(t),0) € H, 5: DC H -> R is given by 6((£,n,<t>)) = D*0(JL) and e = 

°d / e - 

The restrictions C and K of the operators C and K that appear in 
equation (3.1) above to 2/(6), the null space of the operator 6, have 
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natural extensions to bounded operators from V (which is the V-closure 
of N(S)) into V' , the dual of V. The extensions are defined in terms 
of the bilinear forms c(-,-): V x v -♦ R and k(-,-): V x V -» R given by 

A A A A 

(3.5) (C0) (V) - c(<|>,^) - C H (;-<t>00)U-Ni(0) + C d I<D 1 2 <|>,D*^> 0 
and 

(3.6) (K4>)(n 0 - k(i.v) - k H U-$0O)(X-^(JO) + EI<D 2 (J),D*^> 0 

A a 

for $ - (£ , <|>(JD , <t>) € V and y - (X , ^(5t ) , ^) € V. 

The finite element method we develop below could be derived from 
standard energy considerations. While this is not the approach we 
take, it is worth noting that the usual energy expressions can be 
given in terms of the forms, operators and inner products defined 
above. The kinetic energy is given by 

1 

T 0 " ~ <Mu t (t),u t (t)> H , 

the potential or strain energy by 
u 0 ~ — k(u(t) ,u(t)) 

and the Rayleigh dissipation function by 

F 0 “ “ oC u t Ct),u t (t)). 

2 

Written in its weak, variational or distributional form 

<Mu tt (t) ,<|>> H + c(u t (t),4>) + k(u(t),<|>) - <F(t) , <(> > H , 

t > 0, J € V 


(3.7) 
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(3.8) u(0) = 0 u t (0) = 0 

the initial value problem (3.1) - (3.2) in H becomes an initial value 
problem in V'. If we assume that f e L 2 (0,T) and rewrite (3.7), (3.8) 
as an equivalent first order vector system, the theory of abstract 
parabolic systems (see [9], [14]) yields the existence of a unique 
mapping 

u € C( [0,T] ; V) n H x ((0,T) ; V) n C l ([0,T];H) n H*((0,T);V') 

which satisfies (3.7), (3.8). If we are willing to assume further 
that f is Holder continuous then there exists a 

(3.9) u e C( [ 0,T] ; V) n C 1 ((0,T];V) n C l ([0,T];H) n C 2 ((0,T];H) 


A A 

with u(t) + eu t (t) e D, t > 0 which uniquely satisfies the initial 
value problem (3.1) - (3.3). 

In order to demonstrate the convergence of the approximation 
schemes we develop below, we shall require a somewhat more regular 
solution to the initial value problem (3.7), (3.8) than either of the 
conditions on f stated above can guarantee. In addition to (3.9), we 

A 

shall require that u € H 2 ((0,T);V). This can be guaranteed (see [7]) 
if we assume that f e H 1 (-t,T) for some x > 0 with f(t) = 0, t < 0 and 
we modify our original mathematical model so that 

(3.10) F(t) - f(t)0, t € [-t ,T] 


A A 

for some 0 = (O,0(Jl),0), a fixed element in V. We note that with 0 
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chosen appropriately in V, F given by (3.10) may in fact represent an 
improved model of reality when compared with our present choice of F 

/s. 

where 0 = (0,1,0) € H. 

Central to our approach is a cubic spline based Galerkin approxi- 
mation to the initial value problem (3.7), (3.8). For each N =l,2,--- 

N X. 21 

let A denote the uniform mesh (0, — , — , •••, Jt } on [0,51] and let 

N N 

, N N+l 

{B } denote the usual cubic B-splines defined with respect to the 

J J *** 

N N 2 

nodal set A (see [11], [13]). Briefly, each B is a C function on 

J 

Jl X, 

CO, jt] which is a cubic polynomial on each subinterval [(k-1)— ,k— ], 

N N 

N 51 Jl 

k = 1 , 2 , ♦ ♦ • ,N. The support of B is [( j-2)— , ( j+2)— ] n [ 0 , JL ] with 

j N N 

N 51 N 51 N 51 N5LN 

B j (j i ) “ 4 ’ = °' B j ((J±1) i ) = lf and 

^ , r NN+1 N N N N N N 

Defining {0 } by 0 - B -2B -2B_ and 0 = B , J-2,3, • • • ,N+1 , 

J 10 1 1 J J 

N N 

we have 0 (0) - D0 (0) - 0, j - 1,2,---,N+1. With 0 = (1,0,0) and 

i j 0 

~N N N N N N+l 

0 - (0,0 (51), 0 ), j - 1 , 2 , • • • ,N+1 , V = span {0 } is an N+2 

J J J j j=0 

dimensional subspace of V. 

N 

The Galerkin equations in V corresponding to (3.7), (3.8) for 
~N _ N 

u (t) € V are given by 


(3.11) 


(3.12) 


~N ~N .N „N -N „N „N 

<Hu (t),0 > + c(u ( t ) , 0 ) + k(u ( t ) , 0 ) = <F(t) , 0 > , 
tt j H t j j j H 


t > 0, j = 0,1,2, •• • ,N+1 


~N 

u (0) - 0 


.N 

u (0) = 0. 
t 
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Setting 

N+l N .N 

u (t) = I w (t)B , t * 0, 

J J J 

N 

the initial value problem (3.11), (3.12) in V is equivalent to the 
linear , nonhomogeneous , second order N+2 - vector system 

(3.13) 

(3.14) 


d J V7 dw^ 

M® ^-ct) + C N £-(t) ♦ «V(t) - F»(t), 

dt 2 dt 


w (0) = 0 


-( 0 ) = 0 



or, recalling (3.10), by 
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N . N N 

f (t) - f(t)<e,B > - f(t)Cea)B OO + <e,e > ), 

i i H i i 0 


i = 0 , 1 , 2 , • • • ,H+1. 

We consider the sequence of approximating finite dimensional iden- 

— N 

tification problems which consist of finding q € Q which minimizes 

N 

*t 3 2 \i 2 

(3.15) J N (q) - " z(t) l dt 

where for each q € Q, u N (t;q) = (y N (t;q) ,u N (t ,!l;q) ,u N (t, • ;q)) is the 

N 

unique solution to the initial value problem (3.11), (3.12) in V 
corresponding to q = (m T ,E,c D ,m H ,c H ,k H ) 6 Q. In actual practice, for 
a given q € Q, J N (q) is computed as 

j N (q) - It" l w N_i(t ; q) + 4 w ^ t * c i) + w N+i ( ‘ t;CI ' ) " 


N N NT 

where w (*;q) = (w (♦;q),***,w (*;q)) is the unique solution to the 

0 N+l 

N+2 - vector system (3.13), (3.14) corresponding to q € <$. 

With finite dimensional state constraints, the solution of the 
« estimation problem above is, at least in princip le, routine. For 
inverse problems which are closely related to the one we treat here, 
our earlier numerical studies have shown that satisfactory results can 
be obtained using any one of a number of standard computational 
techniques for least squares minimization (for example, Newton's 
method, conjugate gradient, steepest descent, Levenberg-Marquardt , 
etc. , see [2] ) . 
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Our fundamental theoretical result is that each of the approxi- 
mating identification problems and the original problem have 
solutions. Moreover, we show that the solutions to the approximating 
problems, in some sense, approximate solutions to the original 
problem. We require the following lemma. 

LfiUfflta B 1 Suppose {q } C Q with q -» q as N ■* ®. Let u (• ;q ) 
denote the unique solution to the initial value problem (3.11), (3.12) 
corresponding to q and let u(-;q ) denote the unique solution to the 
initial value problem (3.7), (3.8) corresponding to q°. If u(*;q°) € 
H*((0,T);V) then 



as N -> oo . 

Proof 

N 

For each N - 1,2,* •• let P denote the orthogonal projection of H 
N 

onto V defined with respect to the standard inner product on H, 
<*,•>£. Using the approximation theoretic properties of interpolator 
splines, it is not difficult to show that (see [3]) 

(3.17) lim |(P N - I)U,n,<t>) I = 0 

N-»oo I l H 

for each (£,iri,<J>) € H and that 

N A 

(3.18) lim |(P - I)<t> I = 0 

N->® 1 I V 
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for each <J) € V. 

For q = (m^, , E , Cq , nig , Cjj , kg ) € Q it is immediately clear that M, 
c(*,*) and k(*,*), the operator and forms defined in (3.4), (3.5) and 

(3.6) respectively depend upon q. For q° = (m^,E°,cO,m^,c°,k°) € Q 
N N N N N N N 

and q =(m,E,c,m,c,k)€Qwe adopt the shorthand notation 

H° - M(q°), C°( • , • ) - c(q°) ( • , • ) » k° ( • , • ) - k(q°) ( • , • ) , M N - M(q N ), 

c N (*,*) - c(q N )(*,*) and k N (*,*) - k(q N )(*,*). Similarly, we denote 
n *^\T *|T ^ pj 

u(-;q ) and u (*;q ) by u and u respectively. 

From (3.17), the assumption that u° € H l ((0,T);V) and the 
inequality 





it is clear that we need only to consider the first term on the right 
hand side of the above estimate. 

Letting v N (t) - u N (t) - P N u°(t) for t i 0, (3.7), (3.8), (3.11), 
(3.12) and V 1 * c V imply 


(3.19) 


<h n ^,?> h + c n (vN,; n ) + k N ( ; N ,i N ) 


M N rT :n ^ :n 

= <M (I-P )u tt ,$ > H + <(M - M )u tt ,<t> > H 

. j(ff T a o^o :n. js,;,o :n^ 

+ c ((I-P )u t ,4> ) + c (u t ,<|) ) - c (u t ,<t> ) 


+ + - k N cG°,i M ), t . o. ; n € v n 
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(3.20) 

Choosing 


v N (0) - 0 v^(0) = 0. 


i N = v? t (t) € V N , from (3.19) we obtain 


,„n;n 

<H T tf 

’tt’H + 

N v N ) 
t’ v tt ; 




= 

<M N ci-p N )i? t ,^ t . H + 

rw o wNv'o ,;n 
<(M M ) u tt .v tt > H 


+ 

d 

dt 

o 1 'c(i-p h )G°, 

v?)- 

c N ((I. 

-P N )u° t ,^) 


+ 

d 

dt 

- 

c N (u° 

-?)> 

- (c 0 ^ 0 t v N : 

I - C N ^° t ; N )} 

+ 

d 

dt 

k N C(I-P M )G°, 

3>- 

h N ((I 

-P N )U?,V^) 


+ 

d 

dt 

- 


,v?)} 

- <*°^?.v*!) 

- k N C^.vN)} 

- 

d 

dt 

k N c; N ,^) + 

k N (^. 


t > 0. 



Integrating the above expression from 0 to t and recalling (3.20), we 
find 


(3.2i) f <h n ? ,; n > ds + i c N (;;,;;) 

Jq ss’ ss h 2 ' f t' 


l (<m h «-p k )G° s ,^> - 


h 


ss ss 


H 


- c N ((l-P M )u° s ,v^) - Cc°cG2 s .v£) - oFcSS..^)) 

- Aci-p")^) - 

* * N cv*,;»)}ds 


♦ o n ((x-p ,, );°,^) + (o°c;°,;«) - c H <u^)) 
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+ k M (d-p K )i°,v?) + (k°(i°,;?) - * H (i°,;S)) - k N (; 1 ',;;). 


We recall that Q has been assumed to be a closed and bounded subset of 
R® and observe therefore that the forms c°( • , • ) , c N ( • , • ) , k°( • , • ) and 
k N (-,-) are uniformly bounded. These two facts together with the 
repeated application of the inequality 

<a,b> £ |a| |b| £ a|a| 2 + 1 1> | 2 , a > 0 


in (3.21) yield the estimate 


K | ~N | 2 I 

J v ss ds + v: 

J 0 I SSI H I 


,N| 


S Y 0 ( f ( 4F- l«- pH )“ssl! * aft. 


H 


+ 45- c I" , h- , "hI + >1-2.1* + “1-2.1* + 


'H 


‘H 


SS 


+ ra t 7 ^M 2 HoKi 2 )|;° s l!H^ 


2 

V 


* l^u°| v+ |v»| y+ (|^r + |eH-h °|' \%\ 


2 
V 

2 , 2 
V 


2 

V 


^ k <\ * im 


2 

V 
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♦fcPV»l3l v > 



where y q 

is a positive constant. Choosing a > 0 sufficiently small. 

we find 

JqPss (s) I/ s + R (t) l* * a c 

ft f t. 


(3.22) 

)(t) + J 0 ^Cs)ds + yJJ 

where 

o 0 (t) - Y 2 {|(I-P N )u°(t)|^ + | 

(I-P N )u°(t)|^ 



+ k N -q°| 2 C|u°(t) 

+ |u°(t)| ) + |v N (t) 

2 

} 



V 1 t 'v 1 

V 


OjCt) - Tf 3 {|(I-P K )u°(t)|® + 

ci-p N )i° t ct)| 2 




! + > 



and . i = 1,2,3 are positive constants which, do not depend on N. 

a it a vr vr 

Choosing <J> = v^.(t) e V 1 * in (3.19), arguments similar to those 

used above (see [2], [3]) yield 

2 

(3.23) lim |v N (t)| = 0 

N-oo 1 1 v 

for each t € [0,T]. Using u° € H 2 ((0,T);V), (3.18) and an application 
of the Gronwall inequality to (3.22) we obtain the desired result. 


We note that we also obtain 
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(3.24) 


urn Rct)| 2 = 0 

N-*« 1 15 'v 


for each t e [0,T]. From (3.23) and (3.24) we find |u N (t ;q N )-u(t ;q°) I 

>V 

0 and |u^(t;q N )-u t (t;q°) | ■» 0 as H ■* » for each t € [ 0 , T] . 

We remark that it is the L 2 convergence (more precisely, H 
convergence) in (3.16) which necessitates, at least in theory, that we 
be provided with distributed time observations (i.e. observations 
which are continuous in time). It is clear from (3.23) and (3.24) 
that for fits based upon displacement, velocity or slope, time-sampled 
measurements are sufficient. Of course when the approximating 
optimization problems are solved, the integral least squares 
performance indices (3.15) are discretized. Consequently, in 
practice, only discrete measurements of linear acceleration at the tip 
are required. 


Theorem 5,1 Each of the approximating identification problems has a 
solution q^. The sequence {q 1 *} c Q admits a convergent subsequence 
- N 1 - N 1 

(q J } with q -»q € Q as j -»«. if for each q € Q, u( • ;q), the unique 
solution to the initial value problem (3.7), (3.8) corresponding to q, 
is an element in H*((0,T);V) then q is a solution to the original 
identification problem. In addition, the limit point of any convergent 
subsequence of {q N } is a solution to the original identification 
problem as well. 


Proof Standard continuous dependence results for linear ordinary 
differential equations, the fact that Q has been assumed to be a 
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6 N 

closed and bounded subset of R and the form of J are sufficient to 

conclude that a solution q e Q to the N approximating identifica- 
tion problem exists. Once again since Q is a closed and bounded (and 

Q IT 

therefore compact) subset of R , the sequence {q"> C Q admits a 

_N, _N _N, 

convergent subsequence. If {q J } C {q } with q J -* q £ Q as j ■» ® and 
q is any point in Q, then two applications of Lemma 3.1 (the second 
one with the constant sequence {q}) yield 

N N N 

J(q) - lim J *kq S lim J ^(q) = J(q) 

j-*co 

and the theorem is proved. 

Although Theorem 3.1 above guarantees only subsequential 
convergence, in all test and simulation examples we have considered, 
we in fact observe the convergence of the sequence {q N } itself to the 
optimal parameters q. Also, it is not difficult to verify that with 
only minor modification (see [2]) the approximation scheme reported on 
here (together with the convergence theory outlined in the lemma and 
theorem above) could be applied to inverse problems involving the 
estmation of spatially varying parameters (such as linear mass density 
p, flexural stiffness El, or damping coefficient e D I) which appear in 
the equations (2.1) - (2.4). We note of course that when either El or 
o D I are spatially varying, the Euler-Bernoulli equation and 
corresponding boundary conditions are of a slightly different form 
than those given in (2.1) - (2.4) (see [3]). 
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4. Numerical Results 

We used our scheme to attempt to solve the inverse problem which 

was posed above with data obtained from an experiment on the RPL 

structure. We report on our findings and observations here. 

All computer codes were written in Fortran and run on the IBM 

3081 at the University of Southern California. The approximating 

finite dimensional least-squares minimization problems were solved 

using the IMSL implementation of the Levenberg-Marquardt algorithm 

(routine ZXSSQ) , an iterative Newton's method-steepest descent hybrid 

(see[2l). The second order N+2 - vector systems (3.13), (3.14) were 

solved (integrated) in each iteration (for the evaluation of J and 

its gradient) using Gear's method for stiff systems (IMSL routine 

DGEAR) . The integral least squares performance index was approximated 

by a discrete sum over a uniform mesh on [t,,t x ]. The integral inner 

N N N 

products in the definitions of the matrices M , C and K were 

computed using a composite two point Gauss-Legendre quadrature rule. 

N 

The second time derivative of w , or generalized acceleration, 

d‘w N 

— , was computed using a second order centered difference on the 
generalized displacement , 

, d*w N , w N (t+A) - 2 w N (t) + w N (t-A) 

(4.1) (t) * 5 . 

dt 2 tr 

We found this to be a somewhat more stable method for computing 
acceleration (an unbounded measurement) than was a first order 
centered difference on the generalized velocity, 
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(4.2) 


d l w 


N 


dt ’ 


-(t) 


dw N A 

(t + -) 

dt 2 


dw N A 

(t - -) 

dt 2 


A 


Either of the time differencing formulas (4.1) or (4.2) proved to be 
significantly more stable than using the differential equation (3.13) 


d l w 


N 


directly to compute 


dt’ 


-(t) via an inversion of M N . As to why this 


was so, we can only offer the conjecture that the time differencing 
provided, at least to a certain extent, some filtration of the signal. 

Before turning our attention to the experimental data, we tested 
our scheme with simulated data. "True" values for the unknown 
parameters c Q (actually c D I), m H , c H and k H were chosen and a quintic 
spline-based semi -discrete Galerkin scheme applied to the initial 
value problem (3.7), (3.8) was used to generate data. 

Setting p - .03, m T - .15, El = 80.0, l - 4.0 and 

. . 1.0 0 £ t £ 0.05 

f ) == { 

0.0 0.05 < t £ 5.0, 


the fit was carried out based upon observations of linear acceleration 
at the tip at times t ± = .li, i = 2,3, • • • ,50. We note that this is 
equivalent to taking t 0 = .1, t x = 5.0 and using a standard rectangle 
rule with uniform mesh spacing .1 to discretize the integral appearing 
in the definition of the least squares performance index J N . The 
initial estimates c D I = .0035, m H - .035, and k H - .4 were used to 
start the iterative optimization procedure. In (4.1), A was taken to 
be . 1 . Our results are summarized in Table 4 . 1 below . 
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N 

5 N I 

C D X 

-N 

m H 

5 N 

C H 

k N 

«■(?> 

2 

.037537 

.039471 

. 003428 

. 298626 

2.57X10 -1 

3 

. 066997 

.039485 

.003907 

. 298875 

4.37xl0 -2 

4 

. 005063 

.039777 

. 003997 

. 299455 

5.06xl0 -3 

5 

. 005667 

. 039899 

.003971 

. 299787 

7.66xl0 -4 

6 

. 005049 

.040035 

. 004006 

.300087 

4.63xl0“ 5 

True 

value 

.005000 

. 040000 

. 004000 

.300000 


Initial 

Estimate 

.003500 

.035000 

.003500 

.400000 



Table 4.1 


The experiment which we describe below was carried out for us on 
the RPL structure by Dr. Michel A. Floyd, formerly of the Control and 
Flight Dynamics Division of the Charles Stark Draper Laboratory and 
the Department of Aeronautics and Astronautics, MIT. 

The air bearing table was clamped so that the central hub could 
not rotate. The thruster lines for one of the active appendages was 
set to 300 psi and the thruster was fired for .05 seconds (50 milli- 
seconds). With the appendage initially at rest, the firing of the 
thruster was assumed to have begun at time t - 0. Linear acceleration 
at the tip was observed over the time interval 0 to 5 seconds. With a 
sampling period of .005 seconds (5 milliseconds) a total of 1000 
measurements were recorded. The data is plotted in Figure 4.1 below. 
The scale factor for the accelerometer is 5 volts/g (g = 32 ft/sec 2 ). 
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Figure 4 . 1 

The noticeably higher frequency (* 14 Hz) component of the data 
is a torsional mode of the arm excited by the motion of the thruster 
valve mechanisms and inertial and elastic forces applied to the tip of 
the arm by the nitrogen supply hose. The opening or closing of the 
solenoidal valve in the thruster generates an inertial force which 
acts as a torque on the tip of the arm. Consequently, torsional modes 
are excited. Also, in addition to modifying transverse bending 
characteristics, since the hose is attached to the top of the arm, its 
horizontal motion will tend to generate torques which have a 
"twisting" effect. Although the accelerometer is mounted at the 
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center of the arm (and therefore on a nodal line of the longitudinal 
torsional modes, if we assume vertical symmetry) , as the arm twists, 
the accelerometer picks up a component of the earth's gravitational 
force. Since the first torsional mode has a much higher frequency 
than either of the first two flexible modes (.75 Hz and 7.5 Hz, as 
identified from an FFT of the data) and since it is rapidly damped, we 
neglected its contribution to the accelerometer signal, treating it as 
white noise, and left it unmodeled. A detailed discussion of the 
causes of the excitation of the torsional modes and its effect on the 
transverse bending characteristics of the active appendages can be 
found in [6] . 

The physical characteristics of the structure are as follows. 

The arm is made of aluminum and is 4 feet in length, 6 inches in width 
and .125 inches in height. From this we obtain t - 4.0 ft, p - .027 
slug/ft and I - 4.71 x 10 ® (ft)^. The theoretically predicted value 
for E is 15.84 x 10® lb/(ft)^. The mass of the thruster assembly was 
determined to be m^, - .149 slug. From the calibration table in [6], 
we find that a hose pressure of 300 psi is equivalent to a force of 
.297 lb. We set therefore 

- . , 0.297 lb 0 S t S 0.05 

f(t) - { 

0.0 0.05 < t S 5.0 

To serve as a basis for comparison, we neglected the hose effects 
and structural damping (i.e. we chose c d “ m H = °jj “ “ 0) and used 

the standard Euler-Bernoulli model with the parameters p, E, I and m T 
and input f as specified above to generate the plot of linear 
acceleration at the tip given in Figure 4.2. 
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Figure 4.2 


The plot was obtained by integrating the initial value problem (3.13), 

(3.14) with N = 4 and then using (4.1) to compute the acceleration at 

N 

3 2 u 3 2 u 

the free end. The residuals ( (t.JQ ( t , JL) ) over the time 

8t* 3t * 

interval [0,5] are plotted in Figure 4.3. The sum of the squares of 
the residuals (at intervals of .1 seconds) was found to be 3.03. 
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C.OC 1.00 2.00 3.00 -1.00 


t 

I 

Figure 4 . 3 

i 

Using the data on the interval 3.0 to 5.0 (where the contribution 

i 

i from the torsional modes has been significantly damped) with a 
! sampling period of .1 seconds we used our scheme with N - 4 to obtain 
| optimal estimates for the coefficient of viscosity c D and the hose 
| parameters m^, o^ and k^. In the set of runs we are about to describe 
the values of E and m^, were held fixed at their theoretically 
predicted values. A rough calculation based upon "matching" the first 
two observed natural frequencies of the data with the first two modal 
frequencies of the model was used to obtain a crude initial estimate 
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for the ratio k H /m H . Then, using our scheme to minimize over the 
parameters m H and k H only, we obtained the optimal values shown in 
Table 4.2 below. Integrating the system (3.13), (3.14) over the time 
interval [0,5] with m H and k R set to the values in the table and c D = 
c H = 0 the sum of the squares of the residuals (at intervals of . 1 
seconds) was found to be .73. 


m fi (slug) 

k H (lb/ft) 

.039269 

.339935 


Table 4.2 


Next, holding m H and k H fixed at the values shown in Table 4.2, a 
search on c H was carried out (the initial estimate for c H was taken to j 

be zero and c Q was held fixed at zero). Then using the resulting 
values of m H> c H and k H as initial estimates, a fit over all three 
parameters was performed. The result is shown in Table 4.3. The sum ! 
of the squares of the residuals was found to be .728. j 


m H (slug) 

c H (lb*sec/ft) 

k H (lb/ft) 

.043431 

. 004056 

.351385 


Table 4.3 


Continuing to use the same procedure to generate "start up" values, we 
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eventually used our scheme to search over all four parameters c D , m H> 
c H and k H simultaneously obtaining the values given in Table 4.4 and 
the fit plotted in Figure 4.4. The residuals are plotted in Figure 
4.5. The sum of their squares was computed to be .70. 


c D (lb • sec/ (ft) 2 ) 

m H (slug) 

c H (lb'Sec/ft) 

k H (lb/ft) 

127.40 

.0801 

.007804 

.412977 


Table 4.4 



Figure 4.4 
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Figure 4 . 5 

In designing a controller for the RPL experiment, Floyd in [6] 
used model adjustment to tune a simple, undamped, cantilevered beam 
with tip mass model for the active arms (i.e. the arms with the hoses) 
of the structure. He used the following procedure. The air bearing 
table was locked in a stationary position. With the hose depressur- 
ized, an impulsive force was applied to the beam and linear accelera- 
tion at the tip was measured and recorded. Based upon the physical 
assumption that with the hose depressurized, the presence of the hose 
serves only to add mass to the tip of the arm, the parameter m T was 
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adjusted so that the first mode or frequency of the model agreed with 
the first observed cantilever mode (obtained via an FFT) of the data. 
Then, with the hose pressurized, the same experimental procedure was 
carried out. This time however, the modulus of elasticity E of the 
beam was adjusted to compensate for the variation in stiffness which 
results from the presence of the hose. The adjusted values of the tip 
mass, m T , and modulus of elasticity, E, obtained by Floyd are given in 
Table 4.5 below. 


m T (slug) 

E (lb/ (ft) 2 ) 

.254 

17.31 x 10 8 


Table 4.5 


We integrated the system (3.13), (3.14) using the adjusted values of 
m T and E given in the table (and c D “ “ c n “ k jj = 0) and obtained 

the plot shown in Figure 4.6. The corresponding residuals are 
plotted in Figure 4.7. The sum of the squares of the residuals was 
computed to be 5.1. 

Starting with the same basic model, we used our scheme to 
determine the values of m T and E which minimize the sum of the squares 
of the residuals over the time interval [3.0, 5.0] with a sampling 
period of .1 seconds. Taking the theoretically predicted values of m T 
and E (m T - .149 slug, E = 15.84 x 10 8 lb/ (ft) 2 ) as start up values 
for the optimization routine yielded the results given in Table 4.6. 
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The corresponding fit and residuals are plotted in Figures 4.8 and 4.9 
respectively helow. The sum of the squares of the residuals (over the 
interval [0,5]) was computed to be .73. 



0-00 1 .00 2-00 3-00 4.00 

Time (Seconds) 


Figure 4.6 
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Figure 4.9 

In summary, we have seen that analysis of the RPL experimental 
data can be carried out in several ways with a number of different 
models . Our techniques can be used to provide reasonable fits of the 
data to models with or without hose and/or beam damping. Even if one 
attempts to leave the physics of the hose - beam dynamic interaction 
unmodeled and perform "model adjustment" (by adjusting the values of 
the tip mass m T and beam modulus of elasticity E) , our estimation 
techniques provide a much better fit than that obtained using "modal 
matching" methods common in engineering practice. 

One of the primary objectives of our effort here was to 
demonstrate the efficacy of our scheme and in particular, to assess 
its effectiveness when provided with actual experimental data. While 
we are pleased with the results obtained for the RPL data, we are 
careful to point out that to provide a fair and complete evaluation of 
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the usefulness of our models for the RPL experimental structure, a 
more complete and in-depth study involving extensive experimental work 
and statistical analysis would necessarily be required. 
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